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Abstract 

The aim of the paper is to demonstrate the use of the Galerkin method for some 
kind of Volterra equations, determininistic and stochastic as weh. The paper con- 
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apparently well-known results concerning the Volterra equations under considera- 
tion. In the second one we describe a numerical algorithm used and next present 
some examples of numerical solutions in order to illustrate the pertinent features of 
the technique used in the paper. 
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1 Introductiom 



In the paper we investigate a stochastic version of a hnear Volterra equation 
of the general form 

t 

X(t, x)^ j a{t - t)AX{t, x)dT + Xq{x) + f{t, x), (1) 



where t E M+, x G W'', a e Llg^(R+), A is a hnear operator and / some 
mapping. To fix om attention we shall consider the equation (1) in a separa- 
ble Hilbert space H with a scalar product (■, •), a norm | • | and a complete 
orthonormal system {e„}. The equation (1) creates a big class of equations 
and generalizes heat and wave equations and even linear Navier-Stokes equa- 
tion. We refer to the excellent monograph [13] for a rich survey. That kind 
of Volterra equation has been studied by many authors in connection with 
problems arising in mathematical physics, particularly in viscoelasticity, heat 
conduction in materials with memory, energy balance and termoviscoelastic- 
ity. In order to take into account random fluctuations, we have to consider the 
equation (1) with random external force. 

There are our first considerations concerning numerical treatment of stochastic 
Volterra equations, so we will be grateful for readers' remarks and advices. 

Next, we plan to study the probabilistic features of family of trajectories, take 
into account different noises and develope numerical schemes for cases where 
the analytic form of reselvent is not known. 



2 Resolvent approach 

Assume that {fl,J-, {J-'t)t>o,P) is a probability space with a complete right- 
continuous filtration and W{t), i > 0, is a cyhndrical Wiener process with 
values in the space H. Let us omit, for convenience, the space variable x in 
the equation (1) and introduce the process W{t), t > 0, instead of function /. 
Hence, we arrive at the following stochastic Volterra equation 

t 

X{t) = J a{t- T)AX{T)dT + Xq + W{t). (2) 



In this part of the paper we recall some results concerning solutions to (2). We 
restrict our considerations to paper containing so-called resolvent approach to 
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Volterra equation. The notion resolvent or fundamental solution for Volterra 
equation (1) probably comes from Friedman and Shinbrot [5] who studied 
deterministic Voherra integral equations in Banach space. For recent survey 
we refer again to [13]. In the sequel we shall assume that the equation (1) is 
well-posed, that is, that (1) admits resolvent S{t), t >0. 

As in deterministic case, the mild solution to the stochastic Volterra equation 
(2) is of the form 

t 

X{t) = S{t)Xo + j Sit- T)dW{T) , t>Q, (3) 



where S{t), i > 0, is the resolvent family for the equation (1) determined by 
the operator A and the function a. In order to study solution (3) it is enough 
to consider the stochastic convolution 

t 

Ws{t) := J S{t-T)dW{T) , t>0, (4) 



where the stochastic integral is defined according to particular case under 
consideration. 

Stochastic Volterra equations with resolvent approach have been treated by 
several authors, see e.g. [1], [2], [3], [4], [11] and recently [8] and [9]. In the 
first three papers stochastic Volterra equations are studied in connection with 
viscoelasticity and heat conduction in materials with memory. The paper [1] is 
particularly significant because the authors were the first who have extended 
the well-known semigroup approach, applied to stochastic differential equa- 
tions, to the equation (2). The resolvent approach is a natural way of exten- 
sion the semigroup approach which is well-known from the theory of evolution 
equations. That approach enables to follow some results and schemes obtained 
for semigroups. Unfortunately, some results are not valid in our case because 
resolvent family, S{t), t > 0, does not satisfy semigroup property. 

Clement and DaPrato studied stochastic Volterra equation (2) where A was 
self-adjoint, negative operator in the space H, such that 

Ask = -Atfcefc, iJ,k> 0, k eN . 

They considered stochastic Volterra equation (2) driven by the noise term W 
of the form 

+ 00 

{W{t),h)H = T.ih,ek)HPk{t), heH, (5) 
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where {Pk} was a sequence of real-valued, independent Wiener processes. They 
assumed that the kernel function a is completely positive. The consequence of 
completely positiveness of the function a is that the solution s(-,7), 7 > 0, 
to the following equation 

t 

s{t) + ^ J v{t - T)s{T)dT = 1, t>0, (6) 




is nonnegative and nonincreasing for any 7 > 0. In fact, s{t) e [0,1]. In 
[1] regularity of stochastic convolution (4) is studied and holderianity of the 
corresponding trajectories is proved. 

Hypothesis 1 

(i) A is a self-adjoint negative operator and Aek — —/ik^k- 
(a) a is completely positive. 
(Hi) -Tr(A-i) = Er=i(l^A:) < +00. 

Hypothesis 2 There exists 9 e (0, 1) and Cg > such that, 
for all < T < t we have 

t 

j s\ii, a)da < Cei^^-^\t - T\e , 

T 

T 

J[s{ii,T-a) - s{ii,t - (T)fd(T < Ceii'^\t-T\e 


and 

00 

^/x^^<+oo. 

k=l 

Hypothesis 3 There exists M > such that 

|efc(^)|<M, fceN, eeO, 

|Vefe(^)| < M/x^', k&n, eeo, 

where O is a bounded open subset ofM.'^. 
Clement and DaPrato proved the following results. 

Theorem 1 ( [1], Theorem 2.2) Assume that Hypothesis 1 holds. Then for 
any i > the series 
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is convergent in L'^{Q) to a Gaussian random variable Ws{t) with mean and 
covariance operator Qt determined by 



Theorem 2 ( [1], Proposition 3.3) Under Hypotheses 1 and 2, for every pos- 
itive number a < 9/2, the trajectories of Ws are almost surely a-Holder con- 
tinuous. 

Theorem 3 ( [1], Theorem 4- i) Under Hypotheses 1, 2 and 3, the trajectories 
ofWs are almost surely a-Holder continuous in {t,x) for any a e (0, 1/4). 

In the paper [2], white noise perturbation of an integro-differential equation 
arising in the study of evolution of material with memory is studied. In the 
next paper [4], the authors considered evolutionary integral equations as ap- 
pearing in the theory of linear parabolic viscoelasticity forced by white noise. 
As earlier, they studied the stochastic convolution that provides regular solu- 
tions. Additionally, under suitable assumptions the authors proved that the 
samples are Holder-continuous. In the remaining part of the paper [4], the re- 
sults obtained of that paper were put in a wider perspective by consideration 
of equations with fractional derivatives. 

In the paper [3], the authors first proved that Ws{t) is a Gaussian random 
variable for any t > 0. Next, the transition function Pf, t > 0, associated 
with S{t) was considered. When S{t) is the resolvent operator of a stochastic 
Volterra equation, the convolution Ws{t), t > 0, is not a Markov process. 
This fact has the consequence that Pt, t > 0, is not a semigroup and then it is 
not possible to associate to Pt a Kolmogorov equation. However, the authors 
characterized those transition functions such that Pt ip was differentiable for 
any uniformly continuous and bounded function (p. 

There are some other regularity results concerning stochastic convolution (4) . 
In [11] was studied the case when the equation (2) was driven by a corre- 
lated, spatially homogeneous Wiener process W with values in the space of 
real, tempered distributions 5"(M°'). Let T be the covariance of W{1) and the 
associated spectral measure be pL. We considered existence of the solutions to 
(2) in S'{W^) and derived conditions under which the solutions to (2) were 
function- valued and continuous. In that case, the initial value Xq G S"(]R'^), 
a is a locally integrable function and A is an operator given in the Fourier 
transform form 








(7) 



We introduce the following hypothesis. 
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Hypothesis 4 (1) For any 7 > 0, the equation (6) has exactly one solution 
7) locally integrable and measurable with respect to both variables 7 > 
andt> 0. 

(2) Moreover, for any T >0, supjg[o,r] sup^>o k(^)7)| < +00. 
For some special cases the function s{t; 7) may be found exphcitly For instance 

for a{t) = 1, s{t; 7) = e"^*, t>0, 7 > 0; 

for a{t) = t, s{t;^) = cos(^i), t>0, 7 > 0; (8) 

for a{t) = e-\ s{t; 7) = (1 + 7)"^! + -fe~^^+'^^% t>0, 7 > 0. 

In that case, the resolvent family S{t), t > 0, determined by the operator A 
and the function a is given by the formula (7) and has the form 



where r{t) = T ^s{t,(3{-)), t > 0. The following results for stochastic convo- 
lution are consequences of properties of stochastic integral. 

Theorem 4 ( [11], Theorem 2) Let W he a spatially homogeneous Wiener 
process and S{t),t > 0, the resolvent for the equation (2). If Hypothesis 4 
holds then the stochastic equation 



Si<W{t) = J S{t-a)dW{a), t>0 



is a well-defined S'{M.'^) -valued process. For each t > the random variable 
SicW{t) is generalized, stationary random field onW^ with the spectral measure 



UtidX) 



{s{a,P{X)))'da 



li{dX) . 



Theorem 5 ( [11], Theorem 3) Assum,e that Hypothesis 4 holds. Then the 
process S ic W{t) is function-valued for all t > if and only if 

J I J{s{a,(3{X))f da] i^{dX) < +00, t>0. 

If for some e > and all t > 0, 



r 

J J (ln(l + \X\)f+'{s{a,(3{X))fdafi{dX) < +00, 
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then, for each t > 0, S ic W{t) is a sample continuous random field. 

As we have already written, the Volterra equation (1) creates a big class of 
equations. In particular cases, when the operator A and the function a in the 
stochastic Volterra equation (2) are fixed, there is possible to obtain some 
additional regularity results. This is obvious because in particular cases we 
may use some extra features of solutions to (2). For instance, apparently well- 
known is integro differential equation which interpolates heat and wave equa- 
tions, that is the Volterra equation (1), where A — A, the Laplace operator 
and a{t) — t°'~^/r{a), < a < 2, where F is the gamma function. Recently, 
deterministic version was studied in detail by Fujita [6] and, independently, 
by Schneider and Wyss [14] and stochastic version of that integrodifferential 
equation was treated in [7] and [9]. In this paper we shall demonstrate nu- 
merical results obtained for that equation in the deterministic version and the 
stochastic one, as well. 

In the theory of stochastic Volterra equations we are interested not only in 
the existance, uniqueness and regularity of solutions but in some asymptotics, 
too. The paper [8] is concerned with a limit measure of stochastic Volterra 
equation driven by very general noise in a form of a spatially homogeneous 
Wiener process, with values in the space of tempered distributions S'{W^). 
That paper provides necessary and sufficient conditions for the existence of 
the limit measure and additionally, it gives a form of any limit measure. 

Let us summarize the results cited above. The resolvent operators S{t), t > 0, 
corresponding to the Volterra equation (1) do not form any semigroup. There- 
fore it is not possible to obtain such strong results as in the case of evolution 
equations with semigroup generators. In the case of Volterra equations one 
can not use the fractional method of infinite dimensional stochastic calculus. 
That method, used for demonstration of continuity with respect to t for con- 
volutions with semigroups, enables to obtain only some estimates for the con- 
volutions (4). Moreover, it is possible only in some special cases, see e.g. [10]. 
It is clear that the Volterra equation (1), in particular its stochastic version 
(2), is difficult to study. It results from the fact that equation (1) contains 
a wide class of equations. An essential role is played by the kernel function 
which, in general, is assumed to be a locally intcgrable function. As the func- 
tion a and operator A determine the resolvent S(t), t > 0, the type of the 
function a is very important. That is reason why it is so difficult to obtain 
in a general case the continuity of the convolution (4) and some other theo- 
retical results. For more general convolution, significant in many applications, 
hke W'^{t) :— jQS{t — T)'ijj{T)dW{T), where ip is an appropriate process, it 
becomes even far more difficult. 

Therefore, in many cases a numerical support of theoretical (analytical) con- 
siderations is demanded. We need computations for obtaining estimates in 
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regularity results, choosing some paramenters, choosing ftmction a with re- 
quired properties and for visualization of solutions obtained. Numerical anal- 
ysis is particularly important when we arc not able to obtain analytical results. 
Additionally, numerical schemes are especially useful for studying the asymp- 
totics. 



3 Galerkin method for deterministic Volterra equation 



In this section we construct a scheme for numerical solution of the Volterra 
equation (2) without random part, that is, for 

t 

X{t, x) = J a{t - t)AX{t, x)dT + Xq{x). (9) 



We shall consider the case when A is the Laplace operator. Denoting by 
K{x^ t, s) — a{t — s)-^ we can write (9) in the standard form 

t 

X{x, t) = Xo{x) + J K{x, t, s) X{x, s)ds . (10) 



In Galerkin method one introduces the complete set of orthonormal functions 
J = 1, ■ ■ ■ , oo on the interval [0, t], that is fulfilling conditions 

t 

(4(0, Mt)) = / Ur) Mr) dr = 5ij , (11) 



where (•, •) is the scalar product. The set spans a Hilbert space. The 
approximate solution is then postulated in the form of an expansion of the 
unknown true solution in the subspace determined by n first basis functions 

n 

X^{x,t) = J2cj{x)Mt)- (12) 

Inserting (12) into (10) wc obtain 

t 

Xn{x, t) = Xo{x) + j K{x, t, S) Xn{x, s)ds + Sn{x, t) , (13) 




8 



where the function en{x,t) represents the approximation error. Prom (13) we 
have 



£„(X, t) = fn{x, t) - Xo{x) - j K{x, t, S) fn{x, s)ds 



(14) 



From (14) and (12) it can be written as 

n T, n 

Enix, t) = Y^ Ckix) (pkir) - Xo{x) - j K{x, T, s) Ckix) (j)kis) ds . (15) 



k=l 



k=l 



The coefficient functions Cj{x) are determined by the requirement that the 
error function £„(x, t) has to be orthogonal to the subspace 



{(j)j{t),en{x,t)) = j = l,2, ...,n 



(16) 



Then for j = 1,2, ... ,n the following equations hold 



t t r 

Xo{x)(j)j{T)dT -- 


t 



lk=l 



(j)j{T)dT 



- j j K{x, T, s) J2 Ck{x) (f>k{s)ds 

n n k=l 



(pj{T)dT 



The first integral on the r.h.s, due to (12), is very simple 



^Cfe(x) (f)k{T) 
.k=l 



(17) 



n 'j. 

(t)j{T)dT = Ck{x) j (f)k{r)(t)j{r)dr = Cj{x) . (18) 

k=l n 



The second one we calculate in our particular case, K{x,t,s) — a{T — 
as follows 



k=i "'•^ 



" d^Ckjx) 

r1 dx"" 

'==1 LO 



>j{T)dT 

t [ T 

j j - s)(i)kis)ds 



(19) 



(T)rfr . 



Denoting by 
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z z 

9j{x) = j Xo{x)(j)j{T)dT = Xo{x) j (t)j{T)dT , and 



t V T 

O'jk^ j (t>j{T) j a{T - s)(f)k{s)ds 



(20) 
(21) 



LO 

we arrive at the set of coupled differential equations for the functions Cj{x) 



k=l 



dx'^ 



(22) 



This set can be solved numerically (aproximately) by discretization, on a grid 
— {2^1, 3^2, ■ ■ ■ , Xm}- Applying the difference form for the second derivative 



d'^Ckixi) 1 



[ck{xi-i) - 2ck{xi) + Ck{xi+i)], h^Xi- Xi-i 



dx'^ h? 

one obtains from (22) the following set of hnear equations 



gj{xi)^Cj{xi) + ^ %fe [-Ck{xi-i) + 2ck{xi) - Cfe(xi+i)] , 



(23) 



k=l 



where j = 1,2, ... ,n, i = 1,2, ... ,m. Those equations can be written in the 
matrix form 



Ac^g . 



(24) 



Here, {N — n x m) -dimensional vectors c i g have the following structures 



^ ci{x,) ^ 



Cl{Xm) 
C2{xi) 

C2iXm) 

Cn{xi) 

y (^n{Xm) I 



Co 



\CnJ 



^ 9i{xi) ^ 



9l{Xm) 

92{xi) 

92{x^) 

9n{xi) 
^9n{Xm) J 



Go 



\GnJ 



(25) 
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where by Q, Gi we denoted the consecutive m-dimensional blocks of vectors 
c and g, respectively. Then we can write the matrix A in the block form 



^ [An] . . . [A,^' ^ 



A 



[^nl] • • • [^nn] y 



(26) 



where every block is a tridiagonal matrix. The diagonal blocks have the fol- 
lowing structure 



[A, 





1 



J^O'ii l~l~/j2'fl«i 



v 



































' /j2 da 1 ~l~ ^2 Oil / 



(27) 



and nondiagonal ones 



[A 





/l2 Q-jj 


































2 

/j2 dij 











V 






















f^2 "jj 1^2 "ij f^2 'Hj 



(28) 



The set of linear equations (24) can be solved by standard methods, for in- 
stance the LU decomposition [12]. 



4 Stochastic integral 



The essential part of the mild solution (3) is the stochastic convolution (4). We 
present here a particular case when the resolvent S{t) is known in analytical 
form. Let us focus the attention on stochastic Volterra equation (2) with the 
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function a in the form a{t) = t"~^ /r(a). This is the integrodifferential equation 
[6], [8], [13]. For three particular cases, a — 0,1,2 the analytical form of the 
resolvent S is known: 

oo oo 

{S{t)v){x) := j (t)a{t,x-y)v{y)dy= j (pa{t,y)v{x - y)dy, (29) 

— oo — oo 



where the last form comes from the convolution property. We shall illustrate 
the applicability of our numerical algorithms with two cases of the above a 
functions, the case with a — 1 and with a — 2. Then the function (p^ in (29) 
takes the following form 



1 x'^ 

(f)i{t,x)= . exp(— — ) for a — 1 (30) 



(f)2(^t,x) = l{S{t-x) + S{t + x)) for a = 2. (31) 

We assume the process in the form W{t,x) = Wi{t)W2{x). Then, the al- 
gorithm for an approximate construction of the stochastic integral (4) can 
be built in the following way. Let us introduce a time grid {ti — ir : i — 
0, 1, . . . , /} on [0, T], i.e. t = T/I and next a finite sequence of independent 
random variables {Q}, i = 1,2, . . . , I , with standard normal distribution. The 
approximation for the convolution (4) can be written in the form 



7-1 



/ S{t - s) dW{s, x) = ^ S{t - Si)[W{si+i, x) - W{si, x)] (32) 
^=0 

^^^22Ci J (t>a{'t- Si,x-y)W2{y)dy 



For further specification we choose W2{x) = C Xq{x) = C e~^^/^ (the constant 
C represents a 'strength' of stochastic forces). With this assumption, after 
performing the integral (32) for particular 0^ one obtains 



^ 

J S{t- s)dW{s,x) = (33) 



Cr'^y^Ci J exp(— -) for a = l, 

^^'^mr^ ^^A{i+t + Si)' 



and 
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J S{t-s) dW{s, 



X] 



i-i 



i=0 



-{x-t + Sif -{x + t- Sif . 
exp(^ —) + exp(^ —) 



(34) 

for a — 2. 



These explicite forms were inserted into the numerical code. The sequence 
of independent random variables {Ci}-, i = 1, 2, . . . , /, with standard normal 
distribution was generated using subroutines gasdev and rani from [12]. 



5 Numerical results 



We illustrate the efficiency of the numerical approach on two examples of the 
function a, mentioned earlier. As the initial value of the X{t — 0, x) = Xq(x) 
we chose the Gaussian Xo(x) = 1^. The grid in x variable contained m = 
150 intervals with h = 0.2, covering the interval x € [—15, 15]. The dimension 
of the approximation subspace in the Galerkin method was chosen as n = 8. 
The resulting dimension of the matrix A was then 1208x1208 and calculations 
were performed up to T = 6. 

In fig. 1 we show the errors of the numerical solutions to the deterministic 
equation (9) obtained in cases a — 1 (top), and a — 2 (bottom). In both cases 
the errors are relatively small. 

The top part of the fig. 2 displays the solution of the deterministic Volterra 
equation 1 for the case a = 1 as function of time, t e [0, 6]. The bottom part 
of the fig. 2 shows the example of a single stochastic trajectory (i.e. the sum 
of deterministic solution and stochastic integral) for the same case. 

The fig. 3 present the corresponding solutions for the case a — 2. For stochastic 
convolutions the value of the constant C was chosen to be C = 0.1. 
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Fig. 1. Errors of the numerical solutions to deterministic Volterra equation (9) at 
t = 6. Top: a = 1 case, bottom: a = 2 case. The symbols represent the difference 
between the exact (analytical) solution and the numerical solution at given grid 
points. 
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X(t,x) 




Fig. 2. Numerical solution to Volterra equation with a = 1 for t £ [0, 6]: the deter- 
ministic solution (top) and the stochastic one (bottom). 



16 



X(t,x) 




Fig. 3. Numerical solution to Volterra equation with a = 2 for t £ [0, 6]: the deter- 
ministic solution (top) and the stochastic one (bottom). 
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